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Using a Hartree ensemble approximation, we investigate the dynamics of the <p 4 
model in 1 + 1 dimensions. We find that the fields initially thermalize with a Bose- 
Einstein distribution for the fields. Gradually, however, the distribution changes 
towards classical equipartition. Using suitable initial conditions quantum thermal- 
-^-J ' ization is achieved much faster than the onset of this undesirable equipartition. We 

also show how the numerical efficiency of our method can be significantly improved. 

O: 

1 Inhomogeneous Hartree Dynamics 

s/J . In many areas of high-energy physics, e.g. heavy ion collisions and early 

f^ ' universe physics, a non-perturbative understanding of quantum field dynamics 

is required. Computer simulations could fulfill this need, but as is well-known 

the problem is very difficult. Hence one resorts to approximations ,,-Such as 

classical dynamicsa (for recent work seen) , large n or Hartree (see e.g. a) . Here 

we introduce and apply a different type of Hartree or gaussian approximation 

than used previously in which we use an ensemble of gaussian wavefunctions 

""p^ | to compute Green functions. We can only sketch this method here, a detailed 

presentation will appear elsewhere. 

We use the lattice ip 4 model in 1 + 1 dimensions as a test model, which 

Mh. has the following Heisenberg operator equations, 

CD 

43; £ = 7T, iT = (A-n 2 )<p-\<p 3 , (1) 

^ ' with A the lattice laplacian, /i the bare mass and A the coupling constant. 

K^ , Rather than solving the operator equation (Tfl) in all detail, one may focus 

?-h ' on the Green functions. The Hartree approximation assumes that the density 

matrix used to compute these Green functions is of gaussian form, such that all 
information is contained in the one- and two-point functions. The one-point 
function is the mean field and the connected two-point function is conveniently 
expanded in terms of mode- functions, 

(<p x ) = <p x , <^„> conn = £/*#*• ( 2 ) 
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We have restricted ourselves to pure state gaussian density matrices here. The 
Heisenberg equations now provide self-consistent equations for the mean field 
ip and the mode functions f a , 

<p x = A<p x - [p 2 + X(if 2 x + 3 J2 /* /DK 

a. 

ft = A/« - [p 2 + A(3^ 2 + 3 £ /:/r)]/, Q (3) 

a 

In order to simulate a more general density matrix than the gaussian ones 
required in the Hartree approximation, we average over a suitable ensemble of 
Hartree realizations, by specifying different initial conditions and/or coarsen- 
ing in time. In this way we may compute Green functions with a non-gaussian 
density matrix p = J2iPiP?i as 

(0»0v>cann = Z^K&^^L + V { 2 vf] " (disconn.) (4) 

i 

Note that the individual Green functions labeled by i are still computed with 
pure gaussian states pf , as is appropriate in the Hartree approximation. 
The field operator in the Hartree approximation may be written as 

& = Px+53/x6«+/r4]. ( 5 ) 

a 

with a a and a' a time-independent creation and annihilation operators. This 
suggests that the mode functions represent the (quantum) particles in the 
model. It should be stressed that in general tp x is inhomogeneous in space. 
We also note that the equations (0) can in fact be derived from a hamiltonian. 
Since the equations are also strongly non-linear, this suggests that the system 
will evolve to an equilibrium distribution with equipartition of energy, as in 
classical statistical physics. 

2 Observables 

To assess the viability of the Hartree ensemble approximation, we solve the 
equations (H) starting from a number of initial conditions. With Hartree 
dynamics we expect to go beyond classical dynamics, because the width of the 
(gaussian) wavefunction, represented by the mode functions, should capture 
important quantum effects. Of course we cannot expect to capture everything, 
e.g. tunneling is beyond the scope of the gaussian approximation. 

Similarly to using classical dynamics, we expect that after coarse-graining 
in space and time and averaging over initial conditions, we may compute the 
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Figure 1. Particle number densities as a function of u) at various times. The straight line 
gives a BE distribution with temperature T/m = 1.7. (The lattice volume Lm = 32, 
coupling \/m 2 = 0.083 and inverse lattice spacing 1/am = 8; m = m(T = 0)). 

particle number densities and energies from the (Fourier transform of the) 
connected two-point functions of ip and it, 
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)con„ = {n k + \)uJk. 
. (6) 

The over-bar indicates averaging over some time-interval and initial conditions 
as in eq. fl), oj k is the energy, n k the (number) density of particles with 
momentum k and N the number of lattice sites. 

For weak couplings, such as we will use in our numerical work, the particle 
densities should have a Bosc-Einstcin (BE) distribution, 

n k = l/(e^/ T -l), w k = (m 2 (T) + k 2 )3, (7) 

with T the temperature and m(T) an effective finite temperature mass. 

3 Numerical results 

First we use initial conditions which correspond to fields far out of equilibrium: 
gaussians with mean fields that consist of just a few low momenta modes, 
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Figure 2. Left: Particle number densities as a function of ui at times tm = 100, 1400, 2800 
and 4200 (bottom up), starting from BE type initial conditions for the mean field. The 
straight line through the origin is a BE distribution with temperature T/m = 6.7. The 
right figure uses classical dynamics and shows data at tm = 100 and tm = 4200. (Lm = 23, 
\/m 2 = 0.083 and 1/om = 22). 

The aj are random phases and A is a suitable amplitude. Initially the mode 
functions are plane waves, e lkx /v2u)L, i.e. there are no quantum particles 
and all energy resides in the mean field. 

The results in Fig. |l| show that very fast, im^200, a BE distribution is 
established for particles with low momenta. Slowly this thermalization pro- 
gresses to particles with higher momenta, while the temperature T rs 1.7m 
remains roughly constant. Such a thermalization does not happen when us- 
ing homogeneous Hartree dynamics (ip constant in space). The difference may 
be understood, because particles in our method can scatter off the inhomo- 
geneities in the mean field. 

Next we want to speed-up the thermalization of the high momentum 
modes. Therefore we use different initial conditions in which the energy is 
distributed more realistically over the Fourier modes of the mean field. The 
modes are the same as before but mean fields are drawn from an ensemble 
with a BE-like probability distribution, 

P(<p k , n k ) ex exp[-(e^/ T ° - 1)(tt| + «M/2w fc ] (9) 

Now we find density distributions as shown in Fig. ^ (left). Already after a 
short time, tm fa 100 the particles have acquired a BE distribution up to large 
momenta to jm s=s 12. Note that even with these BE-type initial conditions, 
the fields initially are still out of equilibrium: energy is initially carried by the 
mean field only but is quickly, within a time span of tm rs 200, redistributed 
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Figure 3. Particle number densities as a function of u> at times tm ~ 1700 (straight line), 
and 8000, 15000 and 27000 (increasingly curved lines). The drawn lines are fits with Ansatz 
n = co + ci/w. (Lm = 5.7, A/m 2 = 0.083 and 1/am = 22). 



over the modes. 

Since the Hartree dynamics follows from a hamiltonian, one might expect 
that eventually equipartition sets in. To investigate this, we have run a simu- 
lation for long times, with BE-typc initial conditions. The results for several 
times are shown in Fig. 0. The (curved) lines here are fits to the data with 
Ansatz n — Cq + Ci/ui. For times tm^2000, the data show a BE-distribution, 
which was established within £m.;$200. At later times tm ^> 2000 one recog- 
nizes classical equipartition with temperature T — c\ and cq ~ 0. Since the 
energy is gradually distributed over particles with increasingly larger momen- 
tum range, the temperature slowly decreases. 

To see the difference between our Hartree and purely classical dynamics, 
we have repeated the simulation of Fig. (left), now using the classical e.o.m. 
for ip. As can be seen in Fig. @ (right) the initial BE distribution with T = 
6.7m evolves to equipartition much faster. Already at tm w 4200 we find that 
the distribution is well represented by the classical form n — 3.7m/ 'to (the 
curved line in Fig. (right)). 

Finally we try to improve the efficiency of our method. Since we have 
to solve for N mode functions, on a lattice with N sites, the CPU time 
for one time-step grows oc N 2 . However, mode functions corresponding to 
particles with momenta much larger than T should be irrelevant, since these 
particle densities are exponentially suppressed. This suggests that we can 
discard such modes. This is tested in Fig. H, where we compare simulation 
results using all mode functions with results obtained using only a quarter of 
the mode functions. This corresponds with mode functions with initial plane 




Figure 4. Hartrec dynamics with all mode functions (drawn lines) and with a reduced 
number (32 out of 128) mode functions (dots). (Lm = 5.7, \/m 2 = 0.083 and 1/am = 22). 

wave energy w^l7m. Clearly the results for particles with significant densities 
n are indistinguishable. For particles with energy larger than s=s 17m, there 
are no longer mode functions that can provide the vacuum fluctuations and 
consequently the particle density defined by QSj) drops to —1/2. 

4 Conclusion 

We have demonstrated that, using our Hartree ensemble method, we can 
simulate quantum thermalization in a simple scalar field model in real time. 
Only after times much longer than typical equilibration and damping times, 
the approximate nature of the dynamics shows up in deviations from the BE 
distribution towards classical equipartition. See the contribution-.of Smit in 
these proceedings for an estimate of damping times in our model! We have 
furthermore shown that these simulations can be done using a limited number 
of mode functions: only mode functions for particles with energies below a 
few times the temperature are required. 
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